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ABSTRACT 

Given  a  simple  polygon  with  n  sides  in  the  plane 
and  a  set  of  k  point  "sites"  in  its  interior  or  on  the 
boundary,  compute  the  Voronoi  diagram  of  the  set 
of  sites  using  the  internal  "geodesic"  distance 
inside  the  polygon  as  the  metric.  We  describe  an 
C)((n+/:)log(n+^)logn)  time  algorithm  for  solving 
this  problem  and  sketch  a  faster 
O  {{n  +  k)\og{n  +  k))  algorithm  for  the  case  when 
the  set  of  sites  includes  all  reflex  vertices  of  the 
polygon  in  question. 


On  the  geodesic  Voronoi  diagram  of  point  sites  in  a  simple  polygon. 

Boris  Aronovt 
Courant  Institute  of  Mathematical  Sciences 


1.   Introduction. 

Recently  there  has  been  a  significant  upsurge  of  results  concerning 
geometry  inside  a  simple  polygon,  including  an  improved  triangulation  algo- 
rithm [TV],  efficient  algorithms  for  calculating  the  geodesic  center  [PSR]  and 
the  geodesic  diameter  [S2]  of  a  polygon,  a  number  of  new  shortest-path  and 
visibility-related  algorithms  that  require  linear  amount  of  time  beyond  a  tri- 
anguJation  [GHLST],  and  algorithms  for  link  distance  problems  [S],  [Lea]. 
Some  of  the  work  concentrates  on  internal  distance  analogs  of  fundamental 
problems  for  point  sets  in  the  Euclidean  plane.  For  example,  Toussaint  [T] 
developed  an  algorithm  for  computing  the  "relative  convex  hull"  of  a  set  of 
points,  which  is  the  shortest  cycle  containing  all  given  points  and  contained  in 
a  given  simple  polygon,  and  Suri  [S3]  described  an  all-geodesic- furthest- 
neighbors  algorithm  for  simple  polygons.  Our  present  work  extends  another 
fundamental  problem  in  computational  geometry,  namely  that  of  calculating 
the  (nearest-neighbor)  Voronoi  diagram  of  a  set  of  point  sites  (see,  for  exam- 
ple [PrS]),  to  the  case  of  points  in  a  simple  polygon  under  the  interior  geo- 
desic metric  (Figure  1) . 

In  sections  that  follow  we  present  an  algorithm  which  calculates  the  geo- 
desic Voronoi  diagram  of  a  set  S  of  k  points  ("sites")  in  a  simple  polygon 
with  n  vertices,  using  the  length  of  the  shortest  internal  path  between  two 
points  in  the  polygon  as  the  measure  of  distance.  Our  algorithm  runs  in  time 
0{(ji-\-k)\o^{n-\-k)),  which  is  not  far  from  being  optimal  and  which  is  an 
order  of  magnitude  faster  than  a  previous  algorithm  of  [AA]  (see  below). 
Our  work  can  be  considered  as  another  generalization  of  Voronoi  diagrams 
using  a  new  kind  of  metric.  Recent  generalizations  of  similar  nature  include 
[A],  [AE],  [CD],  [F],  [EM],  [LS],  [LW],  [L2],  [OSYl],  [OSY2]. 

One  notable  property  of  the  geodesic  metric  is  that  calculation  of  the 
shortest  path  between  two  points  is  not  an  easy  operation  and,  without 
preprocessing,  must  take  Cl{n)  time  in  the  worst  case.  Also,  a  single  "geo- 
desic bisector"  between  a  pair  of  sites  can  be  the  concatenation  of  0(n)  dis- 
tinct straight  and  hyperbolic  arcs,  which  may  at  first  sight  suggest  that  the 
worst-case  overall  complexity  of  the  geodesic  Voronoi  diagram  is  (d{nk). 
Fortunately,  this  is  not  the  case,  and  we  show  that  the  total  size  of  the 


t  Work  on  this  paper  was  performed  while  the  author  held  AT&T  Bell  Laboratories 
PhD  Scholarship  at  New  York  University. 


Figure  1. 
Geodesic  Voronoi  diagram  of  three  sites. 

diagram  is  only  &{n  +k),  where  only  0(k)  of  the  diagram  vertices  are  "true" 
(degree  3)  vertices.  Roughly  speaking,  the  diagram  can  be  described  by  con- 
structing the  set  of  segments  that  extend  the  edges  of  the  shortest  path  tree 
[GHLST]  from  each  site  s  ^S,  but  only  within  the  Voronoi  cell  of  s;  the  aug- 
mented diagram  can  then  be  regarded  as  the  union  of  all  these  sets  of  exten- 
sion segments,  truncated  to  their  corresponding  Voronoi  cells  and  separated 
from  each  other  by  Voronoi  edges. 

There  has  been  some  work  done  previously  in  design  of  an  algorithm  for 
constructing  the  geodesic  Voronoi  diagram  of  a  set  of  point  sites  inside  a  sim- 
ple polygon.  The  best  result  of  which  we  are  currently  aware  is  the  algorithm 
of  [AA],  which  runs  in  time  0{n k  +  nloglogn  +  k\ogk).  The  portion  of 
their  algorithm  that  necessitates  quadratic  time  is  the  explicit  construction  of 
the  shortest  path  tree  [GHLST]  for  the  full  polygon  from  each  site.  As  these 
trees  are  easily  seen  to  have  a  total  of  n  k  distinct  edges  in  the  simple  case  of 
a  convex  polygon,  the  (worst-case)  quadratic  bound  follows.  We  have  cir- 
cumvented this  problem  by  never  building  the  full  shortest  path  tree  of  the 
polygon  for  each  of  the  sites,  instead  the  tree  (in  fact,  a  different  but  closely 
related  "extension  segment"  structure)  from  each  site  is  constructed  only  for 
the  part  of  the  polygon  that  can  conceivably  lie  in  the  Voronoi  cell  of  the 
site. 

Our  algorithm  uses  a  familiar  divide-and-conquer  strategy  for  obtaining 
the  geodesic  Voronoi  diagram.  However,  peculiarities  of  the  geodesic  metric 
necessitate  a  somewhat  non-standard  implementation  of  this  strategy.    The 


Geodesic  Voronoi  diagram 


-2 


B.  Aronov 


relatively  standard  merge  step,  for  example,  must  be  preceded  by  a  step  that 
extends  a  recursively  computed  diagram  of  the  subset  of  sites  inside  half  the 
polygon  to  the  full  polygon.  The  extension  phase  is  the  least  conventional 
part  of  our  algorithm  and  requires  sweeping  the  triangulation  of  the  polygon 
by  a  polygonal  scan-line. 

Remark:  If  the  set  S  of  sites  contains  all  reflex  vertices  of  P,  a  simpler 
algorithm  can  be  used  to  compute  the  geodesic  Voronoi  diagram  in  time 
o\{n  +k)\og{n  +/:)),  see  discussion  in  Section  5. 

Possible  applications  of  our  algorithm  include  the  closest  pair  problem, 
the  nearest  post  office  problem  and  all-nearest-neighbors  problem  in  the  con- 
text of  a  polygonal  universe,  such  as  a  (polygonal)  island  with  no  interior 
lakes  or  a  polygonal  factory  floor.  It  is  likely  that  other  planar  point  location 
and  proximity  problems  whose  solutions  employ  Euclidean  Voronoi  diagram 
could  be  generalized  to  questions  about  internal  metric  in  a  simple  polygon 
and  could  take  advantage  of  our  algorithm.  For  example,  one  might  wish  to 
investigate  how  the  analogues  of  a  "Delaunay  triangulation"  and  "minimum 
spanning  tree"  behave  in  the  context  of  the  geodesic  metric.^ 

The  paper  is  organized  as  follows:  Section  2  describes  the  general  struc- 
ture of  our  algorithm.  Section  3  examines  some  of  the  more  important 
geometric  properties  of  the  geodesic  Voronoi  diagram.  Section  4  gives  a 
more  detailed  description  of  our  algorithm.  Section  5  outlines  a  simplified 
version  of  the  algorithm,  which  produces  the  diagram  for  a  set  of  sites  that 
includes  all  reflex  vertices.  Section  6  mentions  some  related  open  problems. 

2.   General  outline  of  the  algorithm. 

Our  approach  to  solving  the  problem  is  the  following:  Let  P  be  a  simple 
polygon  with  n  sides,  and  5  be  a  set  of  k  point  sites  inside  P.  The  geodesic 
Voronoi  diagram  of  S  in  P,  denoted  as  Vorp{S),  is  the  partitioning  of  P  into  k 
cells  Vp{si),  •  •  •  ,  Vp{sk)  such  that 

Vp(5,)  =  {A:6P|V;e{l,  •  •  •  ,k]  dp{x,s.;)^dp{x,Sj)}, 

where  dp  is  the  geodesic  distance  inside  P.  Some  initial  processing  includes 
triangulation  of  P  [TV]  and  locating  each  site  of  S  in  the  resulting  triangles 
[ST].  Then  a  balanced  decomposition  [Ch]  of  the  triangulation  tree  is  com- 
puted (see  also  [GHLST]).  Tlie  algorithm  proceeds  by  cutting  P  by  a  chord 
into  two  parts  Pi  and  P2  of  roughly  the  same  number  of  edges  and  recur- 
sively computing  the  diagrams  VorpJ^S-C)  and  Vorp,{S2),  where  5i  {Si)  '^^  the 
subset  of  S  contained  in  P^  (P2,  respectively).  At  this  point,  the  Voronoi 
diagram  Vorp.{Si)  must  be  extended  to  the  full  polygon  P.  The  two  extended 


'  One  generalization  of  Delaunay  triangulation  (briefly  discussed  in  Section  5)  was  inves- 
tigated by  Lee  and  Lin  [LL],  who  defined  it  directly,  and  not  as  the  dual  of  the  Voronoi  di- 
agram. 
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diagrams  thereby  obtained  are  next  merged  in  a  manner  similar  to  the  usual 
Shamos-Hoey  scan  [SH]  (or  its  modification  due  to  Kirkpatrick  [K]).  The 
main  novel  feature  of  our  algorithm  is  the  diagram  extension  step,  which 
itself  takes  0((n  -\-k)\og(n  +k))  time.  The  extension  is  done  by  propagating, 
say,  Vorp^{Si)  from  the  diagonal  cut  e  separating  Pi  and  P2  through  each  of 
the  triangles  in  P2  in  ^  preorder  traversal  of  the  triangulation  (sub)tree  of 
P2-  As  the  construction  proceeds  deeper  into  P?.  two  types  of  events  take 
place:  (i)  a  Voronoi  cell  of  Vorp(Si)  "dies  out,"  or  (ii)  a  vertex  of  P2  is 
reached,  assigned  to  the  Voronoi  cell  of  Vorp(Si)  containing  it,  and  is  incor- 
porated into  the  shortest  path  tree  within  that  cell.  A  "wave  front"  is  main- 
tained and  updated  at  each  event.  At  (i)  a  Voronoi  cell  disappears  from  the 
wave  front  and  new  predictions  are  made  as  for  the  times  of  disappearance  of 
the  two  newly  adjacent  cells.  At  (ii)  the  wave  front  is  split  in  two,  one  for 
each  of  the  unvisited  triangulation  subtrees  rooted  at  the  current  triangle. 
Progress  is  achieved  by  repeatedly  selecting  the  next  closest  event  to  take 
place  and  making  appropriate  changes  in  the  wavefront.  A  more  detailed 
description  can  be  found  in  the  remainder  of  the  paper.  A  sketch  of  a  simpli- 
fied algorithm  for  the  case  when  S  includes  all  reflex  vertices  of  P  is  given  in 
Section  5. 

3.   Properties  of  the  geodesic  Voronoi  diagram. 

One  is  hard  pressed  to  efficiently  compute  the  geodesic  Voronoi  diagram 
before  having  understood  some  of  its  properties  and  having  determined  its 
worst-case  complexity.  In  this  section,  the  definition  of  the  diagram  is  given, 
followed  by  a  list  of  its  basic  properties  and  a  proof  that  it  has  linear  com- 
plexity. 

3.1.    Definitions  and  structural  properties. 

Let  P  be  a  compact  region  in  the  plane  whose  boundary  dP  is  a  simple 
n-gon  with  set  of  vertices  V  =  {vi,  •  •  •  ,v„}.  hcl  S  =  {si,  ■  •  •  ,Sk}CP  be  a 
set  of  k  point  "sites." 

Definition  3.1.  For  any  two  points  v  and  w  of  P,  let  h  (v,w)  be  the  shor- 
test path  between  v  and  w  entirely  contained  in  P.  The  last  vertex  (or  v  if 
there  is  none)  before  w  on  h{v,w)  is  referred  to  as  the  anchor  of  w  (with 
respect  to  v). 

Note  3.2.  As  P  is  closed  and  bounded  by  a  simple  pol^-gon,  h(v,w) 
exists  and  is  indeed  unique  for  every  pair  (v,w)  of  points  of  P.  Moreover, 
h(v,w)  is  piecewise  linear  with  "breakpoints"  at  vertices  of  dP.  For  a  proof 
see,  for  example,  [LP]. 

Definition  3.3.  For  v,  w  ^  P,  let  the  distance  from  v  to  w,  d{v,w),  be  the 
length  of  h{v,w). 

Note  3.4.  Essentially  by  definition,  d(v,w)  is  a  metric  on  P.  Moreover, 
d{v,w)  is  continuous  in  both  \'  and  w  (with  respect  to  Euclidean  metric  on 
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Definition  3.5.   The  Voronoi  cell  of  a  site  5  6S  is 

Vp{s)  =  {xiP  \\lt^S:d{x,s)<  d{x,t)]. 

Note  3.6.   Trivially,   \JVp{s)^P. 
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The  objective  of  our  algorithm  is  to  obtain  the  decomposition  of  P  by 
Voronoi  cells.  Observe  that  applying  known  techniques  for  planar  point  loca- 
tion to  such  a  decorrjposition  allows  efficient  computation  oi  [s  \  x^Vp{s)] 
for  an  arbitrary  query  point  x^P. 

Definition  3.7.  The  shortest  path  tree  of  P  from  a  site  s,  T{P,s),  is  the 
union  of  the  shortest  paths  from  s  to  vertices  of  P.  It  is  known  that  T{P,s)  is 
indeed  a  planar  tree  rooted  at  s  with  straight-line  edges  and  VUl^}  as  the  set 
of  vertices  (see,  for  example,  [LP]).  Assuming  s  is  not  a  vertex,  T{P,s)  has 
exactly  n  edges  and  n  +  1  vertices  and  each  of  its  edges  is  either  a  side  or  an 
interior  chord  of  P  (Figure  2) . 

Definition  3.8.  The  shortest  path  partition  of  P  from  5  is  the  partition  of 
P  into  maximal  sets  each  containing  points  x  all  having  the  same  anchor  with 
respect  to  j  (thus  all  the  paths  h  {s,x)  pass  through  the  same  sequence  of  ver- 
tices of  P). 

Let  e  be  an  edge  of  T{P,s)  and  let  its  endpoint  further  from  s  be  v.  Let 
r  be  the  open  half-line  collinear  with  e  and  extending  from  v  in  the  direction 
of  increasing  distance  from  s.    U  some  initial  section  of  r  is  contained  in  the 


Figure  2. 
Shortest  path  tree  T(P,s)  of  P  from  s  superimposed  on  P. 
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interior  of  P,  we  will  refer  to  the  maximal  such  initial  section  as  the  extension 
segment  of  e  (or  the  extension  segment  emanating  from  v).  Otherwise  e  has  no 
extension  segment. 

Definition  3.9.  Let  the  collection  of  extension  segments  of  edges  of 
T(P,s)  be  denoted  by  E{P,s).  We  will  often  abuse  the  notation  and  write 
E(P,s)  for  the  union  of  the  extension  segments  of  T{P,s).  Note  that  E{P,s) 
contains  at  most  n  segments  — one  per  vertex  of  P. 

Note  3.10.  It  was  shown  in  [GHLST]  that  the  partition  formed  by  split- 
ting P  along  the  segments  of  T(P,s)UE{P,s)  is  a  triangulation  in  which 
points  of  each  triangle  share  an  anchor.  Such  triangles,  however,  are  not 
necessarily  maximal  sets  with  this  property.  If  fact,  it  is  easy  to  see  that  split- 
ting P  along  the  segments  of  E{P,s)  alone  produces  maximal  sets  and  thus 
corresponds  to  the  shortest  path  partition  of  P.  In  particular,  it  is  a  polygo- 
nal partition  of  P  (see  Figure  3). 

Definition  3.11.  If  v,  w  6  P  and  v  ?t  w,  we  will,  following  [PSR],  define 
the  direction  from  v  to  w,  u{v,w),  as  the  unit  vector  at  v  directed  along  the 
first  segment  of  h{v,w). 

Note  3.12.  u{v,w)  is  a  continuous  function  of  v  and  u'  for  all  pairs  of 
(v,w)  except  when  v  ==  w  or  when  there  is  a  segment  of  E{P,w)  emanating 
from  V.  Moreover,  it  is  easily  verified  that  d{v,w)  is  a  continuously  differen- 
tiate function  of  w  with  the  exception  of  points  noted  above  and,  in  fact, 
—  u(w,v)  is  the  gradient  of  d{v,w)  with  respect  to  >v. 


Figure  3.  Shortest  path  partition  of  P  from  s. 
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Definition  3.13.   The  bisector  of  two  (distinct)  sites  s  and  t  is 
bis,t)  =  {x^P  I  d(s,x)^d{t,x)}, 
i.e.,  the  set  of  points  equidistant  from  s  and  t. 

Definition  3.14.  A  set  g  CP  is  star-shaped  around  x  ^P  (with  respect  to 
the  geodesic  metric)  if  Vj  €2  :h{x,y)  C  Q. 

Definition  3.15.   If  s  and  t  are  two  distinct  sites,  let 

H{s,t)  =  {xiP  \his,x)nb{s,t)=^0}. 
In  other  words,  H{s,t)   is  the  maximal  subset   of  P  —  b{s,t)   star-shaped 
around  s.  Note  that  trivially  s  iH{s,t). 

Lemma  3.16.  {H{s,t),b{s,t),H{t,s)]  is  a  partition  of  P.  Moreover, 
for  a  point  a:€/', 

X  ^H{s,t)  <=>  d{x,s)<d{x,t), 

X  ^b{s,t)  <=>  d{x,s)  -  d{xj),  and 

X  iH{t,s)  ■^  d{x,s)>d(x,t). 

Proof.  It  is  sufficient  to  demonstrate  the  equivalence  of  our  original 
definition  of  H{s,t)  and  the  alternative  characterization  given  above. 

Assume  xiH{s,t),  so  h{s,x)nb{s,t)  =  0.  We  must  show  that 
d{x,s)<d(x,t).  Continuously  parametrize  h{s,x)  by  'n:[0,l]^h{s,x)  such 
that  tt(0)  =  s  and  tt(1)=j:.  Note  that  d{'n(T),t)  -  d(u{i),s)  is  a  continuous 
function  of  t  on  [0,1].  It  starts  off  at 

d{'n{0),t)  -d{'Tr(0),s)  =  d{s,t)  >  0. 
It  never  reaches  zero,  as  that  would  indicate  an  intersection  of  h{s,x)  and 
b{s,t).        Thus        it        is        positive        throughout.         In        particular, 
J(7r(l),r)  —  d{'n{l),s)  >  0,  i.e.,  d{x,t)  >d{x,s),  as  desired. 

Conversely,  suppose  xiH{s,t),  so  there  exists  a  point 
y  ih{x,s)nb{s,t).  Thus  d{x,s)  =  d{x,y)  +  d(y,s)  and  d {s ,y)  ^  d (t ,y) .  But 
then  by  the  triangle  inequality 

d{x,t)  ^d{x,y)  +  d{y,t)  =  d{x,y)  +  d(y,s)  =  d{x,s), 
as  desired,  which  completes  the  proof.   □ 

Corollary  3.17.  By  a  simple  continuity  argument,  b{s,t)  actually 
separates  H{s,t)  from  H{t,s)  in  the  sense  that  any  path  from  a  point  of  the 
former  to  a  point  of  the  latter  must  intersect  b(s,t).  In  particular,  H(s,t) 
could  as  well  have  been  defined  as  the  path-connected  (rather  than  star- 
shaped  around  s)  component  of  P  —  b  {s,t)  containing  s  (Figure  4). 

Lemma  3.18.  H(s,t)  U  b(s,t)  is  star-shaped  around  s. 

Proof.  Suppose  h{s,x)(LH{s,t)\J  b{s,t),  i.e.,  h(s,x)nH(t,s)  i^0.  Let 
y  eh{s,x)nH{t,s).  Then 

dis,x)  =  d{s,y)  +  d(y,x)>dit,y)  +  d(y,x)^d{t,x), 
so  X  ^H{t,s),  as  desired.  □ 

Corollary  3.19.  H{s,t)Ub{s,t)  is  a  path-connected  and  simply  con- 
nected subset  of  P. 
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Figure  4. 

H(s,t),  b(s,t),  and  H(t,s)  together  with  the  shortest  path  partition. 

Solid  lines  represent  b(s,t),  dashcd-extension  segments. 

Corollary  3.20.  By  definition  of  the  Voronoi  cell  and  lemma  3.16, 
Vp(s)=      f]     {H{s,t)ijb(s,t)).     In   particular,    lemma   3.18    implies   that 

lis  -  [s] 

Vp{s)  is  star-shaped  around  s. 

For  an  alternate  proof  of  the  fact  that,  in  general,  the  Voronoi  cell  C  of 
a  (point)  site  s  is  always  star-shaped  around  s  in  the  metric  used  to  define  C, 
see,  for  example,  [PrS]. 

Definition  3.21.  P  and  S  are  in  general  position  if  no  vertex  of  P  is 
equidistant  from  two  distinct  sites,  i.e.,  no  bisector  b{s,t)  contains  a  vertex 
of  P. 

Lemma  3.22.  Under  the  assumption  of  general  position,  b(s,t)  is  a 
smooth  curve  connecting  two  points  on  dP  and  having  no  other  points  in 
common  with  dP.  It  is  the  concatenation  of  0{n)  straight  and  hyperbolic 
arcs,  and  the  points  along  b(s,t)  where  adjacent  pairs  of  these  arcs  meet  are 
exactly  the  intersections  of  b{s,t)  with  segments  of  E(P,s)  or  E(P,t).  More- 
over, the  tangent  to  b{s,t)  at  point  x  bisects  the  angle  between  il(x,s)  and 
u(x,t). 

Proof.  Partition  P  into  polygonal  regions  by  removing  from  it 
E{P,s)UE(P,t).  Let  Q  be  the  closure  of  one  of  the  resulting  regions.  By 
construction  there  are  two  anchor  vertices  s  and  r  such  that,  for  any  point  x  in 
Q,  d{s,x)  =  d{s,s)-i-  \s-x\  and  d{t,x)  =  dXt,t)  +  \t-x\,  where  points  are 
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regarded  as  vectors  and  |  •  •  •  |  denotes  the  Euclidean  norm.  Thus  the  con- 
dition x^b(s,t)  is  equivalent  in  Q  to  d{s,s)  +  \s-x\=  d(t,t)  +  \t  -  x\. 
Therefore,  in  Q,  t>{s,t)  is  in  general  a  portion  of  a  hyperbola.  The  hyper- 
bola in  question  has  s  and  t  as  its  two  foci  and,  by  elementary  analytical 
geometry,  its  tangent  bisects  the  angle  between  the  directions  to  s  and  to  t. 
The  following  degenerate  cases  are  possible: 

(i)  The  hyperbola  degenerates  into  a  straight  line  (in  fact,  the  perpendicular 
bisector  of  s  and  t),  if  d{s,s)  =  d{t,t)  and  s  "^  t.  Note  that  this  is  the  only 
possibility  for  a  portion  of  b{s,t)  to  be  a  straight-line  segment  for  P  and 
S  in  general  position.  In  particular,  the  straight-line  portion  of  b{s^t) 
cannot  be  collinear  with  either  of  the  anchors,  as  this  would  imply  s  —  t. 

(ii)  There  are  no  points  equidistant  from  s  and  /  either  because 
\d{s,s)—d{t,t)\> d{s,t)  or  simply  because  the  hyperbola  (or  straight 
line)  in  question  does  not  intersect  Q. 

(iii)  If  s  =  t  and  d{s,s)  =  d{t,t),  a//  points  of  Q  would  belong  to  b{s,t). 
Then,  however,  the  vertex  s  =  t  would  be  equidistant  from  s  and  t, 
violating  the  general  position  assumption. 

(iv)  The  intersection  of  the  hyperbola  with  Q  consists  of  discrete  points. 
This  situation  occurs  only  when  the  hyperbola  in  question  intersects  the 
boundary  of  Q  but  not  its  interior. 

The  above  facts  imply  that  b{s,t)  can  be  regarded  as  the  union  of  (rela- 
tively closed)  straight  and  hyperbolic  arcs  and  discrete  points.  We  now 
proceed  to  prove  a  series  of  properties  of  b{s,t)  which  will  allow  us  to  con- 
clude that  it  is  indeed  a  smooth  curve  connecting  two  points  of  dP  and  com- 
pletely contained  in  the  interior  of  P  otherwise. 

(1)  Every  path-connected  component  of  b{s,t)  meets  dP,  for  suppose  a  com- 
ponent C  does  not— then  there  is  a  closed  curve/?  around  C  that  does  not 
intersect  either  b{s,t)  or  dP.  Thus  p  is  completely  contained  in,  say, 
H{s,t),  implying  that  b{s,t)UH{t,s)  is  not  connected- a  contradiction. 

(2)  There  are  no  non-trivial  cycles  in  b{s,t)  under  the  assumption  of  general 
position,  that  is,  there  is  no  closed  simple  curve  completely  contained  in 
b(s,t).  Suppose  there  exists  such  a  curve  and  consider  a  minimal  cycle  C, 
i.e.,  one  whose  interior  contains  exactly  one  connected  component  of 
P-b(s,t),  say,  H(t,s).  But  then  C  is  a  loop  in  H(s,t)U b(s,t)  that  can- 
not be  contracted  to  a  point,  contradicting  simple  connectedness  of 
H{s,t)Ub{s,t)  (corollary  3.19). 

(3)  Thus  we  have  proven  that  b  (s,t)  is  a  forest  with  each  tree  attached  to  dP 
in  at  least  one  point.  We  will  proceed  to  prove  that  all  leaves  of  b{s,t) 
are  points  on  dP.  Suppose  it  is  not  the  case.  Let  a;  be  a  leaf  of  b  {s,t)  not 
on  dP.  Then  a  sufficiently  small  open  disk  D  centered  at  x  does  not  meet 
dP  and  intersects  b(s,t)  in  a  simple  curve  connecting  x  to  a  point  on  dD. 
Therefore  D  -  b{s,t)  is  completely  contained  in  either  H{s,t)  or  H{t,s). 
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Without  loss  of  generality  assume  that  it  is  contained  in  H{s,t).  As 
H(t,s)[Jb(s,t)  is  star-shaped  around  t,  h{x,t)  must  be  contained  in 
Hlt,s)Ub{s,t),  and  h{x,t)nD  C{H{t,s)Ub(s,t))nD  =  b(s,t)nD, 
which  is  impossible  in  a  non-degenerate  configuration,  as  it  implies  the 
existence  of  a  straight-line  portion  of  b(s,t)  lying  on  the  line  passing 
through  one  of  the  anchor  points,  which  was  shown  in  (iii)  to  contradict 
the  general  position  assumption. 

(4)  There  are  no  isolated  points  of  b(s,t)  on  dP,  for  suppose  there  were 
such  a  point  x.  Then  there  is  a  neighborhood  D  of  x  such  that  D  —  {x}  is 
completely  contained  in,  say,  H{s,t)  which  implies  that  H(t,s)  Ub{s,t)  is 
not  connected,  contradicting  corollary  3.19. 

(5)  No  non-trivial  subsegment  of  dP  can  be  contained  in  b{s,t),  as  that 
would  imply  that  both  anchors  for  a  point  on  this  segment  would  have  to 
lie  on  the  straight  line  containing  the  segment  which,  according  to  (iii), 
contradicts  the  general  position  assumption. 

(6)  The  facts  proven  so  far  imply  that  b{s,t)  is  a  forest  all  of  whose  leaves 
lie  on  dP  and  which  has  no  non-trivial  intersections  with  dP  and  no  iso- 
lated points  on  it.  However,  corollary  3.17  shows  that  removal  of  b{s,i) 
from  P  creates  exactly  two  connected  components.  Therefore,  the  forest 
in  question  consists  of  a  single  path  that  does  not  meet  dP  other  than  at 
either  of  its  endpoints,  thus  proving  that  b{s,t)  is  a  path  connecting  two 
points  on  dP  and  having  no  other  intersections  with  dP. 

Finally,  the  directions  u{x,s)  and  u{x,t)  vary  continuously  along  b{s,t) 
(because  of  note  3.12  and  our  assumption  on  general  position),  so  the  angle 
bisection  property  of  the  individual  hyperbolic  arcs  also  holds  at  their  joining 
points,  implying  that  the  individual  arcs  are  "glued  together"  smoothly.  □ 

Notice  that  b{s,t)  cannot  be  tangent  to  an  extension  segment  of  E{P,s), 
for  it  would  imply  (by  the  previous  lemma)  that  the  directions  from  the  point 
of  tangency  to  s  and  to  /  coincide  so  that  s  =  t,  contradicting  the  general  posi- 
tion assumptions  (cf.  part  (iii)  of  the  proof  above).  In  particular,  each  seg- 
ment of  E(P,s)  intersects  b{s,t)  in  at  most  one  point. 

Note  3.23.  The  fact  that  the  bisector  of  two  sites  is  a  smooth  curve 
whose  tangent  bisects  the  angle  between  the  directions  to  the  two  sites  is 
quite  general.  Not  only  does  it  (obviously)  hold  for  point  sites  in  the 
Euclidean  metric,  but  it  was  also  shown  in  [US]  to  hold  for  arbitrary  convex 
sites  in  the  Euclidean  plane. 

Corollary  3.24.   Under  the  assumption  of  general  position 
interior (Vpis))  =       D     H{s,i)  =-  {x^P  I  Vr^5-{5}  :  d(x,s)<d(x,t)} 

and  Vp{s)  =  closure  {interior  {V p  (s))) ,  where  both  the  interior  and  the  closure 
are  taken  relative  to  P.  In  particular,  Vp{s)  is  a  planar  region  bounded  by 
sections  of  dP  and  of  bisectors  of  the  form  b{s,t),  for  t  ^  s. 
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Corollary  3.25.  Interior  {Vp{s))  is  star-shaped  with  respect  to  s,  as  each 
H{s,t)  is. 

Corollary  3.26.   In  a  non-degenerate  configuration, 

5  9^r  =>  interior {Vp(s))  D  Vp{t)  =  0,  for 
interior(Vpis))CH{s,t)  and  Vp{t)  CH(t,s)U  b(s,t). 

Note  3.27.  The  previous  corollary,  which  excludes  non-trivial  overlap 
among  Voronoi  cells,  is  the  main  reason  for  restricting  our  discussion  to  con- 
figurations in  general  position.  If  there  are  two  sites  equidistant  from  a  sin- 
gle vertex  of  P,  b{s,t)  may  no  longer  be  a  curve,  and  Voronoi  cells  (cf. 
definition  3.5)  may  overlap  non-trivially  (see  Figure  5). 

Lemma  3.28.  Let  v^P,  then  h{s,v)r]b(s,t)  consists  of  at  most  one 
point  (under  the  assumption  of  general  position). 

Proof.  Suppose  there  are  two  such  points  x  and  y.  Without  loss  of  gen- 
erality assume  that  y  is  the  point  further  from  s.  Since  both  H{s,t)  and 
H(s,t)Ub(s,t)  are  star-shaped  with  respect  to  s  (definition  3.15  and  lemma 
3.18)  and  h(s,y)  passes  through  x,  the  subpath  p  of  h{s,y)  from  x  to  y  must 
be  completely  contained  in  bls,t).  By  choosing,  if  necessary,  a  different  x, 
we  may  assume  /?  is  a  straight-line  segment  and  the  two  anchor  points  stay 
constant  over  p.  In  particular,  b{s,t)  contains  a  straight-line  portion  p, 
which  forces  p  to  be  a  part  of  the  perpendicular  bisector  of  the  two  anchor 
points.  However,  p  is  a  portion  of  h{s,y),  and  hence  is  contained  in  the  line 
passing  through  one  of  the  anchor  points  which  in  turn  forces  the  two  anchor 
points  to  coincide.  However,  anchor  points  cannot  coincide  on  the  bisector 
imder  the  general  position  assumption  by  the  proof  of  lemma  3.22.  Contrad- 
iction. □ 


Figure  5. 

A  configuration  not  in  general  position. 

Points  in  the  shaded  region  are  equidistant  from  j  and  r. 
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Corollary  3.29.   If  e  is  a  segment  of  E{P,s),  then  b{s,t)ne  consists  of  at 
most  one  point. 

Proof.    Apply  lemma  3.28,  noting  that  by  construction  e(Zh{s,y),  y 
being  the  endpoint  of  e  furthest  from  s. 

3.2.   Complexity  of  the  geodesic  Voronoi  diagram. 

Since,  by  corollary  3.26,  Voronoi  cell  interiors  do  not  overlap,  the  Voro- 
noi diagram  Vorp(S)  of  a  set  of  sites  5  in  P  is  completely  described  by  the 
union  of  the  Voronoi  cell  boundaries,  the  boundary  of  a  cell  being  the  set 
difference  of  the  cell  and  its  relative  interior  in  P.  Note  that  by  corollary  3.26 
a  point  on  the  boundary  of  a  Voronoi  cell  necessarily  lies  in  at  least  tv.o  cells 
and  thus  is  closest  to  (at  least)  two  sites  simultaneously.  Consider  the  planar 
map  on  P  induced  by  the  union  of  Voronoi  cell  boundaries  with  boundaries 
taken  relative  to  P.  A  point  in  three  or  more  Voronoi  cells,  or  on  dP  and  in 
tvN'O  or  more  Voronoi  cells,  is  called  a  Voronoi  vertex.  Maximal  simple  curves 
contained  in  the  union  of  the  cell  boundaries  and  not  containing  Voronoi  ver- 
tices are  Voronoi  edges.  A  maximal  subcurve  of  a  Voronoi  edge  that  is  a  con- 
tiguous portion  of  a  single  hyperbola  or  a  straight  line  is  referred  to  as  an 
arc.  Endpoints  of  arcs  which  are  not  Voronoi  vertices  are  called  breakpoints. 
Note  that  a  breakpoint  necessarily  has  degree  two,  as  degree  one  is  excluded 
by  corollary  3.24  and  more  than  three  arcs  sharing  an  endpoint  would  consti- 
tute a  Voronoi  vertex.  First,  let  us  consider  the  graph  defined  by  Voronoi 
edges  and  vertices  and  ignore  the  individual  arcs  and  breakpoints.  By  con- 
struction, it  is  a  planar  graph  with  vertices  of  degree  three  and  above.  By 
Euler's  formula  the  complexity  of  this  map  is  linear  in  the  number  of  its 
faces.  However,  for  s^S,  interior (Vp^s))  is  star-shaped  around  s  and,  in  par- 
ticular, is  connected  and  contains  s,  so  Voronoi  cells  correspond  one-one  to 
faces  of  the  planar  map.  Therefore  there  are  only  k  faces  in  this  map,  proving 
that  there  are  0{k)  Voronoi  vertices  and  edges. 

Recall  that  BP  has  n  edges,  thus  adding  n  vertices  and  edges  to  the  com- 
plexity of  the  above  planar  map  as  soon  as  we  "explode"  into  constituent 
segments  the  portions  of  dP  which  play  the  role  of  the  exterior  boundar>'.  It 
remains  to  estimate  the  total  number  of  breakpoints  on  Voronoi  edges.  We 
claim  that  the  number  of  such  points  is  0{n),  thus  bounding  the  total  com- 
plexity of  the  map  in  question  by  0{n+k),  and  proving  that  the  size  of  the 
desired  Voronoi  diagram  is  linear  in  the  size  of  the  input.  Augment  the 
planar  map  induced  by  the  Voronoi  diagram  by  the  following  edges:  For  each 
J 65,  add  to  the  map  segments  of  E(P,s)  truncated  to  interior {Vp(s)).  Seg- 
ments of  E{P,s)  that  do  not  intersect  the  interior  of  Vp(s)  are  simply  dis- 
carded. Note  that  the  only  segments  that  remain  are  those  emanating  from 
vertices  of  P  in  the  interior  of  Vp{s)  (Figure  6). 

A  Voronoi  edge  that  is  a  portion  of  bisector  b{s,t)  in  the  original  map  has 
breakpoints  precisely  at  the  points  of  its  intersection  either  with  segments  of 
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Figure  6. 

Augmented  geodesic  Voronoi  diagram  of  three  sites  of  Figure  1. 

E(P,s)  coming  from  H{s,t),  or  with  segments  of  E{P,t)  coming  from  H(t,s) 
(lemma  3.22).  These  are  exactly  the  portions  of  segments  of  E{P,s)  and 
E{P,t)  that  are  included  in  the  newly  created  structure.  One  should  note  that 
by  star-shapedness  of  interior {yp{s))  the  intersection  of  the  segment  (if  any) 
of  E{P,s)  emanating  from  v  ^  interior  {Vp{s))  with  interior  {Vp{s))  is  indeed  a 
single  segment  emanating  from  v.  Since  such  segments  cannot  intersect  the 
boundary  of  a  Voronoi  cell  at  more  than  one  point  by  corollary  3.29,  and 
moreover,  each  segment  is  in  one-one  correspondence  with  the  vertex  it 
emanates  from,  we  conclude  that  the  total  number  of  breakpoints  is  bounded 
above  by  the  sum  of  the  number  of  vertices  occurring  in  various  Voronoi  cell 
interiors.  By  corollary  3.26  the  interiors  of  Voronoi  cells  are  disjoint,  so  the 
number  of  breakpoints  is  bounded  above  by  n,  thus  ensuring  that  the  sizes  of 
the  original  map  as  well  as  the  new  (augmented)  map  are  linear  in  n+k.  Note 
that  the  augmented  structure  (we  will  refer  to  it  as  the  augmented  geodesic 
Voronoi  diagram  Vorp{S))  is  indeed  a  planar  map,  because  extension  segments 
are  truncated  to  their  respective  Voronoi  cells  and  thus  are  (openly)  disjoint 
from  the  original  map  and  among  themselves.  Our  algorithm  will  actually 
compute  the  latter  map.  Thus  we  obtain: 

Theorem  3.30.  The  complexity  of  Vorp{S)  is  0{n+k)  where  n  is  the 
number  of  sides  of  P  and  |5  |  =  ^.  Moreover,  the  augmented  diagram  Vor*(5) 
also  has  complexity  0{n+k). 


Geodesic  Voronoi  diagram 


-13 


B .  Aronov 


4.    Finally,  the  algorithm. 

We  begin  this  section  with  an  overview  of  our  algorithm  for  computing 
the  (augmented)  geodesic  Voronoi  diagram  of  a  set  of  sites  S  in  a  simple 
polygon  P.  The  steps  that  require  more  detailed  treatment  are  discussed  at 
greater  length  in  the  following  subsections. 

Algorithm  A. 

Input:  A  simple  polygon  P  of  n  sides  and  a  set  S  of  k  point  sites  in  the  inte- 
rior or  on  the  boundary  of  P.  S  and  P  are  assumed  to  be  in  general  position. 

Output:  The  planar  map  induced  on  P  by  Vor'p{S)  (cf.  Section  3.2). 

(1)  Preprocessing 

(i)  Triangulate    P    in    0{n]og\ogn)    time    [TV]     (or    by    a    simpler 

0{n  log  n)  algorithm  of  [GJPT]). 

(ii)  Determine  for  each  site  in  which  triangle  it  lies  by  constructing  the 
planar  map  induced  by  the  triangulation,  preprocessing  it  for  point- 
location  queries  [ST],  and  performing  a  query  for  each  site  in  total 
0((n  +  k)  log  n)  time. 

(iii)  In  0(n\ogn)  time  compute  a  balanced  decomposition  of  the  triangula- 
tion tree  [Ch]  which  will  allow  the  main  body  of  our  algorithm  to 
recursively  cut  the  polygon  into  two  parts  each  having  at  least  one 
quarter  of  the  number  of  sides  in  constant  time  per  cut  (see  also 
[GHLST]). 

(2)  Main  part  of  the  algorithm  (recursive) 

(i)  K  S  consists  of  only  one  site  s,  Vorp{S)  is  a  single  cell  Vp{s)  =  P; 
Vorp{S)  is  then  computed  by  utilizing  the  linear-time  shortest  path 
partition  construction  of  [GHLST] ;  otherwise 

(ii)  If  P  is  a  triangle,  compute  the  Euclidean  Voronoi  diagram  of  S  in 
0{k  log  k)  time,  truncate  it  to  P  in  linear  time,  recording  the  intersec- 
tions of  the  Voronoi  edges  with  dP,  and  return  the  result;  otherwise 

(iii)  Split  P  into  two  roughly  equal  polygons  Pi  and  P^  by  a  cut  and 
divide  S  into  S^  and  Sr  such  that  Si  C  Pi  and  SrCPj^.  This  can  be 
performed  in  constant  time,  as  a  balanced  decomposition  of  P  has 
been  precomputed  and  the  sites  are  already  associated  with  the  trian- 
gles in  which  they  lie,  thereby  making  the  latter  operation  implicit  so 
that  no  processing  is  required  to  partition  S. 

(iv)       Recursively  compute  Vor'p^  (Si)  and  Vor'p^  (Sr). 

(v)  Extend  Vorp^(S,^)  to  Vorp(Sp)  and  Vorp^{Si)  to  Vor'p^Si)  in  time 
0((n  +  k)  log  {n  +  k))  as  described  in  Section  4.1. 

(vi)  Compute  Vorp{S)  by  merging  Vor'p{Si)  and  V<?r^(5^)  in  time 
0{n  +  k).   The  details  of  this  step  can  be  found  in  Section  4.2. 
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In  part  2  of  the  algorithm,  if  either  S^  or  Sj^  is  empty,  omit  the 
corresponding  recursive  call  and  skip  the  merge  phase. 

As  preprocessing  takes  0{{n+k)\ogn)  time  and  the  recursive  part  is 
called  on  two  subproblems  of  sizes  (an,ki)  and  {(\  —  a)n,k2),  where 
ki  +  k')  =k  and  —  <  a  <  — ,  with  the  solution  obtained  from  partial  solutions, 

if  any,  in  time  0((n  +k)\og(n  +k)),  a  simple  recurrence  inequality  bounds  the 
running  time  of  the  algorithm  by  0{{n  +k)\og{n  +k)\ogn). 

4.1.    Extending  the  diagram. 

In  this  section  we  will  describe  a  procedure  implementing  step  2(v)  of 
Algorithm  A.  It  computes  the  augmented  geodesic  Voronoi  diagram  for  a  set 
of  sites  in  a  simple  polygon  P  given  the  diagram  for  the  same  sites  in  a  sub- 
polygon,  where  the  subpolygon  was  obtained  by  cutting  P  along  a  chord,  as 
in  step  2(iii)  of  Algorithm  A. 

Let  P  be  a  closed  region  bounded  by  a  simple  n-gon.  Assume  that  a  tri- 
angulation  of  P  is  available  and  let  e  be  a  chord  of  the  triangulation.  Denote 
(the  closures  of)  the  two  polygonal  regions  into  which  e  cuts  P  by  Pi  and  P2. 
Suppose  5  is  a  set  of  k  sites  in  P 1 .  We  will  describe  an  algorithm  for  extend- 
ing Vorp^{S)  to  P2  and  thereby  obtaining  Vorp{S).  Extension  of  the  aug- 
mented Voronoi  diagram  proceeds  by  traversing  the  triangulation  tree  of  P2 
top-down  starting  from  the  triangle  adjacent  to  e.  "Visiting"  a  triangle  A 
involves  extending  the  augmented  Voronoi  diagram  to  A  and  preparing  the 
structures  for  further  extension  to  unvisited  triangles  adjacent  to  A.  Con- 
struction of  the  augmented  Voronoi  diagram  in  the  interior  of  A  is  performed 
by  a  sweeping  algorithm.  Before  we  describe  the  algorithmic  details,  some 
properties  of  the  augmented  Voronoi  diagram  of  5  in  P 2  need  to  be  exam- 
ined. In  particular,  we  will  show  that  the  edge  structure  of  Vorp(S)  in  the 
interior  of  P 2  is  a  forest  constructible  by  sweeping. 

Consider  the  Voronoi  diagram  of  5  in  P.  Let  G  be  the  graph  formed  by 
Voronoi  edges  and  their  incident  Voronoi  vertices  (refer  to  Figure  7), 
i.e.,  the  edge  and  vertex  structure  of  Vorp{S).  Let  F  be  the  graph  formed  by 
arcs  and  extension  segments  restricted  to  appropriate  Voronoi  cells  together 
with  their  incident  Voronoi  vertices  and  breakpoints  (cf.  Section  3.2),  i.e., 
the  edge  and  vertex  structure  of  Vorp{S).  Note  that  G  and  F  are  planar.  By 
a  slight  abuse  of  notation  we  identify  a  graph  with  its  embedding  in  P  so  that 
vertices  are  identified  with  points  of  P  and  edges— with  curvilinear  segments 
in  P.  In  particular,  we  will  write  G  nP2  (F  nP2)  to  mean  G  (respectively,  F) 
restricted  to  P2,  i.e.,  G  (respectively,  F)  with  vertices  outside  of  P2  deleted 
and  edges  leaving  P2  truncated  down  to  their  points  of  intersection  with  e. 

Lemma  4.1.  G  nP2  is  a  forest  (i.e.,  a  cycle-free  undirected  graph)  with 
leaves  on  6P2. 
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Figure  7.  Solid  lines  represent  G,  together  with  dashed -F. 

Proof.  Suppose  it  contains  a  non-trivial  cycle.  Recall  that  we  identify  a 
graph  with  its  embedding  in  P.  Therefore,  we  may  speak  of  a  minimal  cycle. 
C— one  which  together  with  its  interior  does  not  properly  contain  another 
cycle.  As  G  cannot  have  vertices  of  degree  one  in  the  interior  of  P  (corollary 
3.24),  the  interior  of  C  contains  no  point  of  G,  and  thus  is  a  single  connected 
component  of  P  —  G,  i.e.,  a  cell  of  the  Voronoi  diagram  (strictly  speaking, 
the  relative  interior  of  a  Voronoi  cell).  Thus  there  is  a  Voronoi  cell  whose 
interior  lies  entirely  in  P2-  However,  by  corollary  3.25,  a  Voronoi  cell  must 
contain  its  owner  site  in  its  interior,  contradicting  the  assumption  5  CPi .   □ 

Corollary  4.2.   F  nP2  is  a  forest  with  leaves  on  dP2- 

Proof.  Since,  in  terms  of  its  embedding  in  P,  F  differs  from  G  only  by 
the  extension  segments  truncated  to  their  respective  Voronoi  cells,  it  is  suffi- 
cient to  show  that  extension  segments  can  never  close  a  cycle  in  G.  It  is 
clearly  enough  to  establish  that  no  collection  of  extension  segments-  con- 
tained in  one  cell  can,  together  with  their  incident  endpoints,  form  a  con- 
nected component  that  meets  G  twice.  Since  the  extension  segments  are 
defined  as  portions  of  shortest  paths  from  the  owner  v  of  the  cell,  any  two 
such  segments  are  necessarily  disjoint,  for  otherwise  either  their  common 
endpoint  would  have  two  shortest  paths  to  v,  or  the  segments  would  overlap, 
which  would  imply  existence  of  more  than  one  extension  segment  emanating 
from  a  single  anchor.  It  is  easily  verified  that  the  only  other  way  in  which 
closures  of  two  extension  segments  can  intersect  is  for  one  of  them  to  end  at 
the  anchor  of  the  other  extension  segment  (refer  to  Figure  8). 
In  this  situation,  it  is  convenient  to  consider  the  extension  segment  endpoints 
coinciding  on  dP  to  be  distinct  vertices  (leaves)  of  F.  Note  that  this  situation 
can  only  occur  in  "degenerate"  configurations  when  either  three  vertices  of 
P,  or  two  vertices  of  P  and  a  site  of  S  are  collinear.    Thus  we  may  consider 


'  From  this  point  on  by  an  "extension  segment"  we  will  mean  extension  segment  trun- 
cated to  the  interior  of  the  appropriate  Voronoi  cell. 
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(a)  (b) 

Figure  8. 
Degenerate  extension  segments  of  E(P,s). 
In  both  cases,  e.  is  the  degenerate  segment. 

extension  segments  effectively  disjoint  from  each  other.  Now,  an  extension 
segment  emanates  from  a  vertex  v  of  P  and  terminates  at  the  next  intersection 
with  dP  or  with  a  Voronoi  edge.  Under  the  assumption  of  general  position,  v 
cannot  appear  on  the  bisector  of  two  sites  and  thus  on  any  Voronoi  edge. 
Therefore,  an  extension  segment  intersects  G  at  most  once,  as  desired.  All 
leaves  of  Fn/'2  lie  on  dP2  since  the  same  is  true  of  leaves  of  G  r\P2  and 
extension  segments  create  leaves  on  dP2  only.   □ 

We  now  proceed  to  analyze  the  structure  of  the  forest  F  DPi  in  order  to 
justify  the  sweeping  algorithm  for  computing  it,  which  is  given  at  the  end  of 
this  subsection. 

Lemma  4.3.  Consider  a  triangle  A  which  was  entered  through  its  side  / 
in  pre-order  traversal  of  the  triangulation  tree  of  P2  from  (the  triangle 
incident  to)  e.  Let  x  i  interior  (A)  be  a  point  on  an  edge  g  of  F  nP2-  Then 
the  line  p  tangent  to  g  at  j:  crosses  /. 

Proof.  The  lemma  is  trivially  true  for  an  extension  segment  g  as  g  fol- 
lows the  shortest  path  from  a  site  and  thus  must  enter  A  through  /—there  can 
be  no  anchors  in  the  interior  of  A . 

Consider  an  aic  g  of  b(s,t)  (see  Figure  9). 
Since  the  shortest  paths  from  x  to  s  and  t  enter  A  through  /  and  there  are  no 
anchors  in  the  interior  of  A,  u(x,s)  and  u(x,t)  must  be  directed  in  such  a  way 
as  for  the  half -lines  containing  them  to  cut  /.  But  by  lemma  3.22  p  bisects  the 
angle  between  u{x,s)  and  u{x,t),  so  it  must  also  meet  /.  □ 

Lemma  4.4.  Let  A  and  /  be  as  above.  Then  all  edges  of  G  nP2  that 
meet  /  intersect  it  transversally. 

Proof.  Let/'  and  /"  be  the  two  remaining  sides  of  A.  Consider  an  arc 
g  of  b(s,t)  non-transversally  meeting  /  at  x.  Note  that  x  is  an  interior  point 
of /since,  on  the  assumption  of  general  position,  b{s,t)  cannot  contain  ver- 
tices of  P.  Since  the  shortest  paths  from  j:  to  jr  and  to  /  must  not  come  from 
/',  /",  or  the  interior  of  A,  both  u{x,s)  and  u(x,t)  must  lie  in  the  closed 
half-plane  bounded  by  /  and  not  containing  A.  The  only  possible  way  of 
keeping  this  consistent  with  lemma  3.22  that  states  that  the  tangent  to  g  at  x, 
i.e.,  /,  must  bisect  the  angle  between  ii{x,s)  and  u{x,t)  is  for  the  latter  two 
vectors  to  coincide  and  point  along  /.  But  coincidence  of  u(x,s)  and  u{x,t)  on 
a  bisector  is  impossible  in  general  position  by  the  proof  of  lemma  3.22. 
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Figure  9. 

Contradiction.  □ 

Note  4.5.  The  obvious  extension  of  the  previous  lemma  to  all  edges  of 
F  nP2  does  not  necessarily  hold.  The  reason  is  that,  though  ordinarily  exten- 
sion segments  meet  /  transversally,  in  some  degenerate  configurations  (cf. 
Figure  8)  an  extension  segment  might  coincide  with  interior  (f).  We  will 
hereafter  refer  to  these  as  degenerate  extension  segments. 

We  now  proceed  to  define  a  "natural"  orientation  on  FPiPt?  formaliz- 
ing the  intuitive  notion  "away  from  e",  that  transforms  Fr\P2  into  a  root- 
directed  forest.  Orient  edges  of  F  nP2  ^  follows:  Again,  let  A  be  a  triangle 
of  P2  entered  through  /.  Given  an  edge  g,  a  point  j:  on  ^  in  the  interior  of  A 
or  of  /,  the  tangent  p  to  g  at  x  intersects  /  transversally.  Consider  the  direc- 
tion of  traversal  of  g  for  which  the  velocity  vector  at  x  points  along  p  away 
from/Hp  (Figure  10a).  In  case  x  C/ direct  the  velocity  vector  inward  A  (Fig- 
ure 10b). 

A  degenerate  extension  segment  that  overlaps  /  is  oriented  so  as  to  emanate 
from  its  anchor  point.  We  claim  that  this  orientation  is  well-defined  and 
indeed  makes  FHPi  a^  root-directed  forest.  To  show  this,  it  is  sufficient  to 
prove  that  the  above  method  indeed  assigns  a  unique  orientation  to  each  edge 
independent  of  the  point  x  selected  to  define  the  orientation  and  that  no  ver- 
tex has  out-degree  greater  than  one. 

Lemma  4.6.  The  orientation  of  an  edge  g  is  well-defined,  i.e.,  there  is  a 
traversal  of  g  consistent  with  the  above  requirements  on  the  direction  of  the 
velocity  vector  at  every  point  of  g. 

Proof.  It  is  clear  that  the  orientation  of  a  degenerate  extension  segment 
cannot  be  inconsistent  as  it  is  defined  globally  rather  than  locally.  Consider 
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Figure  10.  Definition  of  orientation. 

any  other  type  of  edge  g.  Let  us  start  by  restricting  our  attention  to  a  triangle 
A  entered  through  its  side  /.  The  orientation  is  defined  in  terms  of  where  the 
tangent  to  g  intersects  /.  As  the  tangent  varies  continuously  along  g  (lemma 
3.22)  and  continues  to  intersect  /  (lemma  4.3),  it  continues  to  meet  f  on  the 
same  side  of  the  point  of  tangency,  so  the  orientation  of  g  is  consistent  on  a 
single  connected  portion  of  g  in  the  interior  of  A .  The  ambiguous  case  of  the 
tangent  coinciding  with  /  is  excluded  by  lemma  4.4.  Thus  it  is  sufficient  to 
show  that  the  orientation  stays  consistent  when  crossing  chords  separating  tri- 
angles. And  indeed  it  is,  for  suppose  g  crosses  the  edge/'  separating  A  from 
A',  i.e.,  A'  is  entered  though  /'.  By  lemma  4.4  g  intersects  /'  transversally. 
By  smoothness  of  bisectors  (lemma  3.22),  the  tangent  to  g  at  g  Hf  still 
meets  /.  According  to  our  definition  of  orientation,  g  will  be  oriented  out- 
ward from  A  in  the  interior  of  A  near  /',  it  will  be  oriented  from  A  to  A'  at 
its  intersection  point  with  /'  and  it  will  be  oriented  inward  A'  in  the  interior 
of  A'  near/',  thus  making  the  orientation  consistent  everywhere  along  g.  □ 

So  the  orientation  of  an  edge  is  well-defined.  Let  us  proceed  to  classify 
the  vertices  of  the  forest  FnP2  and  to  prove  that  the  orientation  indeed 
yields  a  root-directed  forest.  We  observe  that  the  leaves  of  the  directed 
forest,  i.e.,  vertices  with  zero  in-degree  are  exactly  the  points  of  intersection 
of  F  and  e  together  with  vertices  of  3P2  that  have  extension  segments  of  F 
emanating  from  them.  These  vertices  of  F  nP2  indeed  have  degree  one  and 
their  incident  edges  are  directed  away  from  them.  An  internal  vertex  of 
F  nP2  (i.e.,  neither  a  root  nor  a  leaf)  in  the  interior  of  A  has  out-degree  of 
exactly  one,  as  shown  by  the  following  lemma.  The  case  of  an  internal  ver- 
tex landing  on  a  chord  of  the  triangulation  will  be  handled  separately  (cf. 
note  4.9  below). 

Lemma  4.7.  Every  vertex  v  of  F  0^2  in  the  interior  of  A  has  out- 
degree  one. 
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Proof.  Suppose  there  are  two  or  more  edges  emanating  from  v.  Note 
that  these  edges  are  Voronoi  edges  and  thus  v  is  necessarily  a  Voronoi  vertex, 
for  an  extension  segment  cannot  emanate  from  a  point  in  the  interior  of  P. 
So  suppose  there  are  two  or  more  Voronoi  edges  emanating  from  v.  Choose 
a  pair  of  adjacent  edges.  As  the  edges  are  adjacent,  there  must  be  a  single 
Voronoi  cell  Vp{t)  (locally)  "wedged"  between  them  and  the  said  edges  must 
be  portions  of  bisectors  b{s,t)  and  b(t,r),  for  some  sites  s  and  r.  Let  t  be  the 
anchor  of  v  with  respect  to  t.  Note  that,  by  definition  of  orientation,  the 
tangents  /^  to  b(s,t)  and  I2  to  b{t,r)  at  v  intersect  /when  extended  back 
through  V.  By  construction,  Vp{t)  is  (locally)  wedged  between  the  half-lines 
of  /i  and  I2  extending  away  from /(see  Figure  11a). 

In  particular,  by  star-shapedness  of  V/>(/),  r  must  lie  in  the  wedge  W,  contrad- 
icting the  fact  that  all  anchors  are  contained  in  the  polygon  over  which  the 
diagram  has  already  been  constructed  and  thus  must  lie  on  /  or  on  the  side  of 
/  opposite  A  and  be  visible  from  v  through  /—in  particular,  /  must  lie 
between  v  and  t. 

Suppose  there  is  no  edge  emanating  from  v.  If  v  is  a  breakpoint  (i.e.,  a 
meeting  point  of  a  bisector  and  an  extension  segment),  two  arcs  of  the  bisec- 
tor join  smoothly  at  v,  thus  guaranteeing  the  presence  of  at  least  one  outgo- 
ing edge.  Hence  v  must  be  a  Voronoi  vertex.  As  G  nP2  has  no  leaves  in  the 
interior  of  A,  there  must  be  two  or  more  Voronoi  edges  incident  to  v.  Con- 
sider the  tangents  to  these  edges  at  v  and  choose  the  two  edges  whose 
tangents  intersect  /  closest  to  either  of  its  two  endpoints  (cf.  Figure  lib). 
Similarly  to  the  above  argument,  the  larger  wedge  at  v  bounded  by  the  t\vo 
edges  must  (locally)  belong  to  a  single  cell  Vp{t)  for  some  t  iS,  and  the  said 
edges  must  be  portions  of  b{s,t)  and  b{t,r),  respectively,  for  some  s,r  iS. 
Let  /i  (/2)  be  the  line  tangent  to  b{s,t)  {b{t,r),  respectively)  at  v.  Let  W  be 
the  smaller  wedge  at  v  bounded  by  the  half-lines  of  li  and  1 2  that  do  not 
intersect/.  By  star-shapedness  of  V/>(r),  u(y,t)  must  point  into  W,  again  con- 
tradicting the  fact  that  the  path  h(v,t)  exits  A  through  /  and  its  portion 


b(t.r) 


Figure  11. 
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contained  in  /is  a  straight  line.   □ 

The  vertices  of  F  DP 2  on  dPi  that  are  neither  anchor  points  nor  lie  on  e 
are  roots  as  shown  by  the  following  lemma. 

Lemma  4.8.  Let  A  be  a  triangle  entered  through  /  and  let  its  remaining 
two  sides  be  /'  and  /" .  No  vertex  v  of  F  n  P2  in  the  interior  of  /'  or  /' '  has 
an  outgoing  edge  g  that  emanates  from  v  into  the  interior  of  A. 

Proof.  The  edge  g  in  question  must  intersect  (say)  /'  transversally  at  v 
(otherwise  it  is  a  degenerate  extension  segment,  and  so  coincides  with  /'  by 
lemma  4.4  and  note  4.5)  and,  by  lemma  4.3  and  smoothness  of  arcs,  the 
tangent  to  ^  at  v  must  intersect  /,  making  g  an  incoming  edge  at  v  by  defini- 
tion of  orientation  on  g.   □ 

Note  4.9.  There  can  exist  internal  vertices  of  F  HP 2  that  lie  on  a  chord 
/'  separating  two  triangles  A  and  A'  of  P2  with  A  entered  through  /and  A' 
through  /'.  In  this  case,  lemma  4.8  shows  that  there  can  be  no  edge  leaving 
V  into  A  and  reasoning  similar  to  lemma  4.7  shows  that  there  is  exactly  one 
edge  emanating  from  v  into  A',  thereby  showing  that  v  has  indeed  out-degree 
one  and  thus  is  a  valid  internal  vertex  of  the  root-directed  forest. 

Corollary  4.10.  F  nP2  is  a.  root-directed  forest  with  leaves  on  e  and  at 
anchor  points,  roots  at  non- anchor  vertices  of  F  on  dP2  —  e,  all  internal  non- 
root  nodes  being  Voronoi  vertices  and  breakpoints  and  having  in-degree  of  at 
least  two. 

Definition  4.11.  Let  A  be  a  triangle  entered  through  /  and  let  /'  and  /" 
be  its  two  remaining  sides.  Define  the  sweep -curve  s^  to  be  the  broken  line 
consisting  of  a  segment  parallel  to  /  at  distance  t  from  it  extending  from  /'  to 
/"  together  with  two  portions  of  /'  and  /"  from  the  endpoints  of  /to  the 
endpoints  of  the  said  segment.  Let  r  be  the  distance  from /to  the  vertex  of  A 
incident  to  /'  and  /" . 

Corollary  4.12.  In  this  notation,  s^  (for  0<T<r)  separates  the  vertices 
of  FnP2  irito  two  sets  in  such  a  way  that  all  (directed)  edges  between  two 
vertices  of  different  sets  emanate  from  below  s^  (i.e.,  from  vertices  lying  on 
the  side  of  s^  that  contains  f)  and  terminate  above  s^.  Moreover,  if  t  is  such 
that  s^  passes  through  a  vertex  v  of  F  nP2  in  the  interior  of  A  or  /,  there  is 
precisely  one  edge  incident  from  v  and  the  said  edge  extends  from  v  into  the 
area  above  s^  (Figures  12a&b). 

Lemma  4.13.  FDA  (  =(F  Pi Fo)  ^^ )  can  be  constructed  by  sweeping  A 
with  s^. 

Proof.  All  that  is  needed  to  build  F  in  A  is  to  systematically  locate  a 
pair  of  children  with  a  common  parent  in  F  PI  A,  and  replace  the  children  by 
the  parent,  thereby  advancing  the  sweep-curve  through  the  forest  one  step. 
This  in  turn  reduces  to  maintaining  the  intersection  of  F  with  current  sweep- 
curve  and  repeatedly  locating  the  pair  of  adjacent  edges  whose  intersection 
point  lies  on  ^^  with  least  t.   □ 
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Figure  12.  Behavior  at  a  vertex  of  FC\P^. 

We  are  now  ready  to  describe  the  algorithm  for  extending  'Vor''p,{S)  to 
Pi-  Throughout  the  algorithm  F  r\P^^,rreni  for  the  current  polygon  P current  is 
maintained  while  P current  grows  from  Pi  to  P.  More  precisely,  maintain  the 
planar  map  induced  by  dPcurrem  together  with  FHPcurrent  (i-C-,  ^o^'pc^^rA^))- 
For  each  bounded  face  R  of  this  map  (we  will  refer  to  it  as  a  region),  main- 
tain the  following  information:  the  site  s  to  whose  Voronoi  cell  R  belongs, 
the  anchor  s  for  the  shortest  path  from  5  to  an  arbitrary  point  of  R,  and  the 
value  d(s,s).  Therefore,  given  an  arbitrary  face  R  of  the  planar  map  and  a 
point  X  in  R,  the  site  s  closest  to  x  and  the  distance  d{s,x)  can  be  determined 
in  constant  time. 

We  will  assume  that  the  planar  map  for  Pi  is  available  initially.  Since 
possessing  Vorp^{S)  implies  knowledge  of  the  sorted  order  in  which  edges  of 
F  meet  the  edge  e  separating  P2  from  Pi,  we  can  build  in  linear  time  a 
search  tree  containing  this  information.  More  precisely,  the  tree  will 
represent  those  regions  of  the  induced  planar  map  on  Pj  that  are  adjacent  to 
e.  Namely,  it  will  contain  for  each  such  region  R  the  owner  of  R  (the  site  to 
whose  Voronoi  cell  R  belongs),  the  anchor  of  R  (the  last  vertex  on  the  shor- 
test path  from  the  owner  to  an  arbitrary  point  in  R)  and  the  weight  of  the 
anchor  (the  distance  between  the  owner  and  the  anchor).  The  tree  must  sup- 
port insert,  delete  and  split  operations  in  time  logarithmic  in  its  maximum  size 
(see,  for  example,  red-black  trees  of  [GS]). 

The  algorithm  for  constructing  Fr)P2  from  Vorp^{S)  (and  thus  FflPi) 
and  thereby  obtaining  F  will  proceed  by  performing  some  initial  setup  at  the 
edge  e  and  then  traversing  the  triangulation  tree  of  Pn  starting  from  the  tri- 
angle adjacent  to  e: 

Algorithm  B 

Input:  A  triangulated  polygon  P  separated  by  a  chord  e  of  the  triangulation 
into  subpolygons  Pj  and  P2,  and  the  planar  map  induced  on  Pi  by  Vor'p,  (5). 
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Output:  The  planar  map  induced  on  P  by  F  and  dP  (i.e.,  Vor*p{S)). 

(1)  Initial  step. 

(i)         From  Vorp^{S),  extract  the  ordered  list  L  of  regions  adjacent  to  e. 

(ii)        Convert  the  list  into  a  search  tree  T. 

(iii)  For  each  pair  of  adjacent  edges  bounding  a  single  region  on  L  (and 
thus  in  T)  compute  their  point  of  intersection  without  checking 
whether  the  resulting  intersections  are  feasible  in  any  sense.  Perfoim 
a  point  location  query  on  each  of  the  points  obtained  and  collect  them 
into  buckets  corresponding  to  the  triangle  where  each  of  those  inter- 
section points  lands  (an  equally  feasible  approach  would  be  to  priori- 
tize the  buckets  immediately). 

(2)  Process  every  triangle  of  the  triangulation  of  P2  in  ^  pre-order  traversal 
of  its  triangulation  tree,  starting  at  the  triangle  adjacent  to  e  and  treating 
the  tree  as  rooted  at  that  triangle.  The  order  in  which  children  are 
visited  is  irrelevant.  Upon  entering,  through  edge/,  a  triangle  A,  whose 
remaining  edges  are  /'  and  /",  whose  vertex  shared  by  /'  and  /"  is  v, 
and  whose  associated  tree  T  represents  the  regions  of  the  planar  map 
intersecting/,  perform: 

(i)  Let  X  be  the  common  endpoint  of  /  and  /' .  Locate  the  anchor  y  of  the 
region  to  which  x  belongs.  Then  consider  the  ray  out  of  x  directed 
away  from  y.  If  this  ray  partially  overlaps  /,  create  a  degenerate 
extension  segment  that  emanates  from  x  and  follows  /  until  its  first 
intersection  with  another  edge  of  F  as  represented  by  T.  If  the  said 
ray  intersects  the  interior  of  A,  create  an  extension  segment  emanat- 
ing from  X  and  directed  away  from  y.  In  either  case,  add  the  intersec- 
tion of  the  new  extension  segment  with  its  neighbor  in  T  to  the 
bucket  of  the  appropriate  triangle  and  create  a  new  region  sandwiched 
between  the  new  extension  segment  and  /'  with  its  owner  being  the 
owner  of  y,  its  anchor  being  x  and  the  weight  of  its  anchor  being  the 
sum  of  d{x,y)  and  the  weight  of  y.  Otherwise  there  is  no  need  to 
create  a  new  extension  segment  or  a  new  region.  Repeat  the  same 
procedure  for  the  other  endpoint  of  /.  Note  that  creation  of  a  region 
involves  insertion  into  T. 

(ii)  Construct  a  priority  queue  Q  containing  all  intersections  in  the  bucket 
associated  with  A.  g  is  ordered  by  increasing  distance  from  the  line 
containing  /,  which  amounts  to  sweeping  across  A  with  s.^.  (Let  us 
note  that  a  different  ordering  would  work  just  as  well  as  long  as  on 
every  straight  segment  in  A  intersecting  /  the  relation  reduces  to  the 
total  order  along  the  segment. 

(iii)  The  intersection  p  in  Q  closest  to  /is  then  repeatedly  deleted  from  Q. 
If  p  is  actually  valid,  i.e.,  it  refers  to  an  intersection  of  two  curves 
currently  adjacent  on  s.  (i.e.,  bounding  a  common  region  in  T),  it 
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must  correspond  to  the  disappearance  of  the  region  between  them 
from  T  and  the  replacement  of  the  two  bisector  arcs  (or  a  bisector  arc 
and  an  extension  segment)  by  a  new  bisector  arc.  At  this  point  the 
region  is  deleted  from  T.  If  p  does  not  represent  a  valid  intersection, 
it  is  simply  discarded.  The  two  newly  created  pairs  of  curves  adjacent 
on  Jt  have  their  intersections  computed  and  located  in  the  triangula- 
tion.  If  they  land  in  A,  they  are  added  to  Q;  otherwise,  they  are 
added  to  the  appropriate  buckets.  The  process  repeats  until  the 
queue  is  empty. 

(iv)  When  the  construction  inside  A  is  completed,  the  two  trees 
corresponding  to  /'  and  /"  are  produced  by  splitting  the  current  ver- 
sion of  T,  (conceptually)  cutting  s^  at  v.  Note  that  s-  coincides  with 
f'Uf"  at  this  point,  so  splitting  T  correctly  creates  the  two  initial 
trees  representing  sq  for  the  two  unvisited  triangles  adjacent  to  A.  In 
case  /'  (respectively  /")  is  actually  an  edge  of  P,  the  corresponding 
search  tree  is  converted  to  adjacency  information  for  F  along  it. 

We  proceed  to  analyze  the  time  complexity  of  Algorithm  B.  Steps  l(i) 
and  l(ii)  clearly  take  linear  time.  As  there  is  only  a  linear  number  of  initial 
candidate  pairs  of  adjacent  curves  on  s-,  step  l(iii)  requires  0{{n+k)\ogn) 
time.  One  can  easily  see  that  the  time  complexity  of  part  2  is  0{\ogn)  per 
candidate  intersection  for  point  location,  0{\og(n+k))  for  queue  mainte- 
nance operations  (as  no  candidate  is  ever  in  two  queues),  0(\og{n+k))  for 
tree  deletion  for  each  valid  intersection  (i.e.,  a  Voronoi  vertex  or  a  break- 
point), and  0(\og(n+k))  for  tree  insertion  and  split  per  vertex  of  P2-  The 
last  three  boimds  follow  as  soon  as  it  is  shown  that  the  size  of  all  queues  and 
trees  is  at  all  times  0{n+k).  In  fact,  each  candidate  intersection  is  either 
created  initially  (there  are  only  0{\L  \)  —  0{n+k)  of  such)  or  added  during  a 
deletion  from  or  insertion  into  a  tree  (at  most  two  new  candidates  per  inser- 
tion, one  new  candidate  per  deletion).  There  are  at  most  n  insertions  and 
0{n+k)  deletions— note  that  the  universe  does  not  necessarily  decrease  all 
the  time;  in  fact  it  can  grow  by  as  much  as  n  due  to  introduction  of  new 
extension  segments,  which  can  occur  at  most  once  per  vertex  of  ^2-  Thus  the 
total  size  of  all  queues  and  trees  is  at  every  point  of  the  algorithm  bounded 
by  0{n+k),  as  desired.  In  other  words,  part  2  spends  logarithmic  time  per 
candidate  intersection  which,  there  being  0{n+k)  candidate  intersections, 
bounds  the  execution  time  of  part  2  by  0{{n  +k)\og{n  +k)). 

Therefore  the  extension  of  the  augmented  Voronoi  diagram  of  a  set  of  k 
points  from  a  portion  of  an  n-gon  to  the  whole  polygon  is  accomplished  in 
O {{n  +  k)\og{n  +  k))  time  by  Algorithm  B,  providing  the  desired  bound  for 
part  2(v)  of  Algorithm  A.  Observe  that  it  is  impossible  to  improve  Algo- 
rithm B  so  that  it  runs  in  time  0{n-^k)  without  modifying  the  remaining 
parts  of  Algorithm  A,  for  such  an  improvement  would  reduce  the  time  com- 
plexity of  the  entire  algorithm  to  O {{n  +  k)\ogn) ,  contradicting  the  lower 
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bound  discussed  in  Section  6. 

4.2.   Merging  the  two  diagrams. 

In  this  section  we  describe  a  linear-time  algorithm  for  merging  the  Voro- 
noi  diagrams  of  two  sets  of  sites  in  a  polygon  assuming  that  a  chord  of  the 
polygon  separates  the  two  sets  of  sites.  Again,  before  detailing  the  algorithm, 
we  will  examine  some  properties  of  Voronoi  diagrams  of  such  sets  of  sites. 

Definition  4.14.  Let  A  be  a  set  of  sites  in  P  and  j:  be  a  point  in  P.  Then 
the  distance  from  j:  to  A  is  defined  by  d{x,A)  =  min  d{x,s).  Let  a  shortest 

path  from  x  to  A  be  h{x,A)  =  h{x,s)  where  j  €A  is  such  that 
d(x,s)  =  d{x,A).  Note  that  h{x,A)  is  not  uniquely  defined  for  ;c  on  an  edge 
of  Vorp (A). 

Definition  4.15.  Let  A,  B  CP  be  disjoint  non-empty  sets  of  sites.  The 
bisector  of  A  and  5  is  b{A,B)  =  {x  €  P  \d{x,A)  =  d(x,B)}. 

Note  4.16.  It  is  easily  checked  that  b(A,B)  is  the  union  of  the  Voronoi 
edges  of  Vorp  (A  UB)  that  separate  the  cell  of  an  A-site  from  that  of  a  B-site 
and  the  Voronoi  vertices  that  lie  on  the  boundary  of  at  least  one  cell  of  a  site 
in  A  and  at  least  one  cell  of  a  site  in  B.  In  particular,  b{A,B)  has  complexity 
linear  in  the  sum  of  sizes  of  A,  B,  and  P  since  it  is  (the  union  of)  a  collection 
of  Voronoi  edges  and  vertices.  Moreover,  for  any  Voronoi  edge  e  that 
b(A,B)  contains,  it  also  contains  the  two  Voronoi  vertices  incident  to  e. 
Hence  b{A,B)  is  a  subgraph  of  G  (as  defined  in  Section  4.1).  Moreover, 
b{A,B)  cannot  terminate  at  Voronoi  vertices  in  the  interior  of  P  or  contain 
isolated  Voronoi  vertices,  as  the  fact  that  a  Voronoi  vertex  v  belongs  to 
b(A,B)  already  implies  that  v  is  in  Voronoi  cells  of  at  least  one  A-site  and  at 
least  one  5-site.  This  in  turn  implies  the  existence  of  at  least  two  (and  at  least 
one  for  v  6  dP)  Voronoi  edges  emanating  from  v  that  separate  cells  of  sites 
from  different  sets.  In  fact,  this  argument  shows  that  there  is  always  an  even 
number  of  edges  belonging  to  b(A,B)  and  emanating  from  an  arbitrary  Voro- 
noi vertex  in  the  interior  of  P.  Thus  b(A,B)  can  be  decomposed  into  edge- 
disjoint  simple  paths  and  cycles.  Paths  connect  two  points  of  dP  and  cycles 
can  be  considered  not  to  touch  dP,  for  a  cycle  that  does  touch  it  can  be 
regarded  as  a  path  connecting  a  point  of  dP  to  itself. 

Definition  4.17.  //(A,5)  =  {x^P  \d{x,A)<d{x,B)}.  In  other  words, 
H(A,B)  is  the  set  of  points  in  P  lying  closer  to  some  site  of  A  than  to  any  site 
ofB. 

Notice  that  trivially  {H{A,B),  b{A,B),H(B,A)}  is  a  partition  of  P  and, 
by  a  simple  continuity  argument,  b(A,B)  actually  separates  H(A,B)  from 
H{B,A)  in  the  sense  that  they  cannot  be  connected  by  a  path  that  does  not 
meet  Z?(A,B). 

Lemma  4.18.  Suppose  s  ^A.  Let  Vp(s)  denote  the  Voronoi  cell  of  s  in 
Vorp  {A)  and  V'p{s)  denote  its  cell  in  Vorp  (A  UB).   Then 
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V'p(s)  =  Vpis)r\{H(A,B)Ub{A,B)). 


Proof. 


V'p(s)  =  {x\d(s,x)^    min    dir,x)} 

r  iA  UB 

—  {x  \d{s,x)  =  min  d{r,x)]  fl  [x  \  min  d{r,x)  <  min  d{r,x)] 

r  i.  A  r  (.  A  r  d  B 

=  Vp{s)f^{x\d{x,A)^d{x,B)]  =  yp(5)n(//(A,5)U^(A,fi)).n 

Similarly  one  shows  that  interior{V' p{s))  =  interior {yp{s))r\H{A,B).  In 
particular,  since  interior {V p{s))  is  connected,  it  is  simply  the  connected  com- 
ponent of  the  set  interior {Vp{s))  —  b{A,B)  containing  s. 

Lemma  4.19.  Suppose  that  Vorp(A)  and  Vorp{B)  are  given.  Then, 
given  a  point  on  every  connected  component  of  b{A,B),  both  b{A,B)  and 
Vor'p(AUB)  can  be  constructed  in  time  C>(|A  |+ |B  |-+-/i)  where  n  is  the 
number  of  vertices  of  P. 

Proof.  From  each  point  given  on  b{A,B)  build  connected  components  of 
b{A,B)  by  the  usual  Shamos-Hoey  scan  of  the  two  Voronoi  diagrams  (see, 
for  example,  [L]  or  [SH]).  Note  that  one  can  follow  a  bisector  of  two  sites  as 
long  as  the  two  anchor  points  stay  fixed.  The  anchor  points  change  precisely 
when  the  bisector  crosses  an  extension  segment,  thereby  moving  from  one 
region  of  the  shortest  path  partition  around  a  site  to  another.  At  this  point 
one  can  recover  the  new  anchor  point  as  it  is  simply  the  anchor  associated 
with  the  shortest  path  partition  region  being  entered.  The  fine  structure  of 
the  Voronoi  cells  (namely  their  shortest  path  partition)  is  somewhat  similar 
to  the  "spokes"  of  Kirkpatrick  [K].  A  complication  arises  when  vertices  of 
degree  greater  than  two  are  encountered  on  b(A,B).  At  such  a  Voronoi  ver- 
tex V  the  list  of  owners  and  anchors  of  adjacent  Voronoi  cells  (and  shortest 
path  partition  regions)  is  known,  allowing  one  to  determine  all  the  edges 
leaving  v  that  separate  a  cell  of  an  A-site  from  that  of  a  fi-site,  thus  enabling 
the  Shamos-Hoey  scan  to  trace  all  branches  of  b{A,B).  The  additional  effort 
required  is  proportional  to  degree  of  v  in  Vorp{A[JB),  thus  it  is  linear  over 
the  construction  of  b(A,B).  Note  that  this  complication  could  be  avoided 
completely  by  taking  "general  position"  to  mean,  in  addition,  that  no  four 
points  are  on  a  geodesic  circle,  thus  excluding  Voronoi  vertices  of  degree 
above  three  and  forcing  degree  of  at  most  two  for  b{A,B). 

Now  split  both  augmented  Voronoi  diagrams  along  b{A,B)  and  collect 
the  portions  of  Voronoi  cells  reachable  from  their  owners;  the  resulting  col- 
lection is  (by  lemma  4.18)  the  augmented  Voronoi  diagram  of  AUB.  Note 
that  one  need  not  be  concerned  with  updating  the  shortest  path  partition  of 
the  cells,  as  the  effect  of  the  merge  on  it  is  limited  to  truncation  of  extension 
segments  down  to  new  cell  boundaries,  which  is  done  automatically  during 
the  Shamos-Hoey  scan.  □ 
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Observe  that  the  two  steps  described  above  need  not  be  done  separately. 
It  is  more  convenient  to  perform  truncation  of  Voronoi  cells  by  b  (A, 5)  at  the 
time  it  is  being  constructed,  concurrently  updating  the  two  Voronoi  diagrams 
and  merging  them  along  b(A,B). 

Note  4.20.  To  repeatedly  locate  the  intersection  of  a  (simple)  arc  of  a 
bisector  with  the  boundary  of  "current"  regions  of  Vorp{Si)  and  Vorp{S2)  in 
an  efficient  manner,  we  follow  Kirkpatrick  [K]  and  introduce  a  further 
refinement  of  our  map  that  will  guarantee  that  each  individual  region  has 
bounded  complexity,  thereby  allowing  us  to  trace  an  arc  of  the  bisector 
within  each  such  subregion  in  constant  time.  In  this  extra  refinement,  in  addi- 
tion to  the  extension  segments,  every  vertex  of  dP  and  every  Voronoi  vertex 
lying  in  a  Voronoi  cell  are  connected  by  a  shortest  path  to  the  owner  of  the 
cell  (actually,  one  needs  only  connect  each  such  point  to  its  anchor  by  a 
straight  segment,  so  the  resulting  structure  is  still  linear).  The  structure 
described  coincides  with  the  "Voronoi  triangulation"  of  [AA]  (Figure  13). 
Moreover,  since  each  segment  introduced  is  contained  in  the  shortest  path  to 
the  owner,  it  will  not  intersect  the  "contour"  b(A,B)  more  than  once.  It  is 
also  easy  to  see  that  every  face  of  the  resulting  map  is  bounded  by  two  seg- 
ments and  (in  general)  an  arc  of  a  hyperbola  or  a  polygon  boundary  segment, 
and  thus  has  complexity  bounded  by  a  constant.  It  is  easy  to  verify  that  this 
finer  structure  allows  Kirkpatrick's  tracing  procedure  to  run  in  linear  time. 


Figure  13. 

Voronoi  triangulation  of  P  for  three  sites  of  Figure  1. 
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Corollary  4.21.  If  all  components  of  b(A,B)  meet  dP,  Vorp{A  UB)  can 
be  produced  from  Vorp{A)  and  Vorp(B)  in  time  0(\A  |-f  \B  \  +  n). 

Proof.  It  is  sufficient  to  locate  all  points  of  b(A,B)ndP.  Since  both 
Vorp{A)  and  Vor'p{B)  are  linear  in  size,  it  is  possible  to  trace  dP,  in  linear 
time  partitioning  it  into  0(n+k)  straight  segments  each  of  which  lies  in 
exactly  one  region  of  each  partition.  Each  segment  contains  points  not  only 
nearest  to  a  unique  s  ^A  and  a  unique  t  C  B  but  also  having  a  unique  pair  of 
anchors  (s,t)  with  d(s,s)  and  d{t,t)  known.  Thus,  on  each  segment,  the  ana- 
lytic expressions  for 

d{x,A)  =  d{x,s)  —  \x  —  s\  +  d(s,s)  and 
dix,B)  =  d{x,t)^  \x-t\+dit,t) 
are  known.  Since  the  two  functions  are  well-behaved  (one  would  like  to  say 
"polynomial  of  degree  two",  but  they  in  fact  involve  radicals),  the  points 
where  they  coincide  are  bounded  in  number  by  a  constant  (actually,  by  two) 
and  can  be  computed  in  constant  time.  Each  such  point  corresponds  to  an 
intersection  of  b{A,B)  and  dP,  thereby  enabling  us  to  compute  all  such  inter- 
sections in  linear  time  and  permitting  the  application  of  lemma  4.19. 

Note  that  the  case  when  the  functions  d{x,A)  and  d(x,B)  coincide  is 
easily  shown  to  violate  the  general  position  assumption  as  in  such  a  case  s  =  t. 

D 

Lemma  4.22.  If  e  is  an  internal  chord  of  P  and  A  and  B  are  sets  of  sites 
of  P  on  opposite  sides  of  e,  b{A,B)  contains  no  cycles.  In  particular,  the  con- 
ditions of  corollary  4.21  are  satisfied. 

Proof.    Denote  the  part  of  P  that  contains  A  hy  Pi  and  the  complemen- 
tary part  by  P2  (refer  to  Figure  14). 
Notice  that  H{A,B)  is  the  union  of  interiors  of  Voronoi  cells  of  A-sites  (in 


Figure  14 
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Vorp(A  DB))  together  with  the  Voronoi  edges  separating  these  cells.  Hence, 
if  there  is  a  cycle  C  in  b(A,B),  it  must  completely  enclose  a  cell  of,  say,  A- 
site  s,  i.e.,  there  is  no  path  in  H(A,B)  connecting  s  to  dP.  Choose  x  to  be  the 
point  of  e  closest  to  s.  Observe  that  either  x  lies  in  the  interior  of  e  and 
u(x,s)  is  perpendicular  to  e,  or  x  is  one  of  the  endpoints  of  e  and  the  angle 
between  u{x,s)  and  the  vector  pointing  along  e  away  from  x  has  measure 
>  tt/2.  Extend  h{x,s)  past  s  until  it  intersects  dP  at  point  y  (Figure  15). 
We  claim  that  h(s,y)CH{A,B),  contrary  to  our  assumption  that  s  is  enclosed 
by  a  cycle  CCb{A,B).  Consider  z^h(s,y)  and  t  ^BCP2-  If  x^h(z,t), 
d(z,t)  =>  d{z,x)  => d{z,s)  (note  that  s  ^  t,  since  A  and  B  are  disjoint),  so 
let  us  assume  x(.h{z,t).  By  the  above  observation,  the  angle  between  u(x,z) 
{—u{x,s))  and  u{x,t)  has  measure  >:Tr/2.  It  was  shown  by  [PSR]  that  the 
side  of  a  geodesic  triangle  lying  against  such  an  angle  is  strictly  the  longest. 
In  particular, 

d{z,t)>d{z,x)>d{z,s). 

Hence  d{z,B)>  d(z,s),  and  therefore  z  ^H{A,B),  as  asserted.   □ 

Corollary  4.23.  The  diagrams  Vor*p{Si)  and  Vorp{Sji)  in  step  2(vi)  of 
Algorithm  A  can  be  merged  in  time  0{n+k). 

To  summarize,  the  procedure  for  merging  the  two  diagrams  consists  of 
locating  all  points  of  b{A,B)  on  BP  by  a  linear  scan  followed  by  tracing  of 
every  connected  component  of  b(A,B)  a  la  Kirkpatrick  [K]. 

5.   A  special  case. 

It  is  interesting  to  discuss  a  variation  of  our  algorithm  that  computes  the 
geodesic  Voronoi  diagram  for  a  set  of  k  sites  among  which  are  all  reflex  ver- 
tices   of   the   enclosing    Aj-gon    in    O  ((n  +  k)\og(n  +  k))    time.     This    variant 
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exploits  the  following  property  of  the  geodesic  Voronoi  diagram  of  such  a  set 
of  sites— the  Voronoi  cells  are  star-shaped  in  the  Euclidean  metric,  as  a  bend 
on  a  path  from  a  point  in  a  cell  to  the  owner  of  the  cell  would  produce  a 
reflex  vertex,  i.e.,  a  closer  site.  The  general  strategy  is  still  divide-and- 
conquer,  with  the  polygon  recursively  split  into  two  roughly  equal  parts  and 
the  resulting  diagrams  merged  a  la  Shamos  and  Hoey  [SH].  However,  the 
extension  step  is  avoided,  rather  a  structure  similar  to  the  extension  of  Voro- 
noi diagram  into  "site-less"  regions  is  built  at  the  bottom  of  recursion  and 
maintained  together  with  the  Voronoi  diagram  proper  through  the  merging 
steps.  The  rationale  behind  this  approach  is  that  the  only  sites  that  can  own  a 
given  point  of  the  polygon  are  those  visible  from  it.  Consequently,  during 
the  merge  step,  the  only  points  of  one  half  of  the  polygon  that  can  be  "taken 
over"  by  sites  of  the  other  half  are  those  visible  from  the  sites  through  the 
"window"  of  the  dividing  chord.  In  particular,  it  is  enough  to  construct  the 
Voronoi  diagram  of  the  set  of  sites  of  one  half  of  the  polygon  as  visible 
through  this  window.  This  notion  is  identical  to  the  notion  of  the  "peeper's 
Voronoi  diagram"  of  [BS],  but  fortunately  in  our  case  it  is  easily  seen  to  have 
linear  size.  Moreover,  we  can  show  that  the  geodesic  Voronoi  diagram  of  a 
set  of  sites  in  a  polygon  with  "peeper's"  Voronoi  diagrams  attached  to  some 
of  the  polygon  edges  (in  fact,  to  those  which  are  chords  of  triangulation  of 
the  original  polygon)  is  a  structure  linear  in  the  size  of  the  polygon  and  the 
number  of  sites  and  two  such  structures  can  be  merged  in  linear  time.  This 
yields  an  O  {{n  +  k)\og{n  +  k))  algorithm  for  constructing  the  geodesic  Voro- 
noi diagram  in  this  special  case.  The  details  of  this  faster  algorithm  will  be 
reported  separately. 

Note  that  the  algorithm  sketched  above  allows,  in  particular,  to  compute 
the  geodesic  Voronoi  diagram  of  the  vertices  of  a  simple  n-gon  in  time 
Oi/iXogn).  An  alternative  0{n\ogn)  algorithm  is  provided  by  work  of  Lee 
and  Lin  [LL]  (see  also  [Chew]  and  [WS]),  which  involves  the  calculation  of 
the  generalized  Delaunay  triangulation  of  a  simple  polygon  P.  It  is  a  triangula- 
tion of  the  polygon  with  the  property  that  the  circle  circumscribed  around 
each  face  of  it  does  not  contain  (in  its  interior)  any  vertex  of  P  visible  simul- 
taneously from  all  three  vertices  of  the  said  triangular  face.  In  the  case  of  a 
convex  polygon,  this  definition  gives  precisely  the  conventional  Delaunay  tri- 
angulation of  P,  which  is  the  dual  of  the  (conventional)  Voronoi  diagram  of 
the  vertices  of  P,  which  in  turn  coincides  (inside  of  P)  with  the  geodesic 
Voronoi  diagram  of  the  vertices  of  P.  In  fact,  it  is  not  difficult  to  observe 
that,  for  a  general  simple  polygon  P,  the  generalized  Delaunay  triangulation 
is  (essentially)  the  dual  of  the  geodesic  Voronoi  diagram  (in  fact,  it  can  be 
shown  to  have  more  edges  than  the  dual).  [LL]  (and  also  [Chew],  [WS])  pro- 
vide an  0{n\ogn)  algorithm  for  constructing  this  triangulation,  and  the 
diagram  can  be  reconstructed  from  it  in  linear  time  by  visiting  the  neighbors 
of  each  vertex  in  cyclical  order  and  thereby  reconstructing  its  Voronoi  cell. 
Thus  we  obtain  an  alternative  0{n\ogn)  algorithm  for  computing  the  geodesic 
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Voronoi  diagram  of  a  set  of  vertices  of  a  simple  polygon  (with  respect  to  that 
polygon). 

6.   Open  problems. 

It  remains  an  open  problem  to  determine  whether  the  geodesic  Voronoi 
diagram  of  a  general  set  of  k  sites  in  an  n-gon  can  be  computed  in  time 
O  {{n  +  k)\og{n  +  k))  or,  for  that  matter,  any  faster  than 
n((/i  +k)]og{n  +k)\ogn).  Let  us  observe  that  the  Euclidean  closest-pair  prob- 
lem is  linear-time  reducible  to  the  geodesic  Voronoi  diagram  problem  as  fol- 
lows: Given  a  set  5  of  it  points,  construct  in  linear  time  a  rectangular  box  B 
large  enough  to  enclose  the  convex  hull  of  S.  At  this  point,  computing 
VorgiS)  produces  the  (conventional)  Voronoi  diagram  of  S  truncated  to  B. 
However,  it  has  been  shown  (see,  for  example,  [PrS])  that  the  shortest 
Euclidean  distance  between  two  points  of  S  is  realized  by  a  pair  of  points 
whose  cells  are  adjacent  and  the  segment  connecting  which  intersects  the  edge 
common  to  the  two  cells— hence  the  intersection  point  must  lie  in  the  convex 
hull  of  S  and  thus  in  B.  In  particular,  we  can  locate  the  closest  pair  in  linear 
time  given  the  truncated  Voronoi  diagram,  which  shows  an  Cl(k\ogk)  lower 
bound  for  computing  the  geodesic  Voronoi  diagram.  Recalling  the  @{n-\-k) 
bound  on  the  complexity  of  the  geodesic  Voronoi  diagram,  we  obtain  the 
best  currently  known  lower  bound  of  n{n  ■+k\ogk). 

There  are  a  number  of  other  open  problems  quite  closely  related  to  geo- 
desic Voronoi  diagrams.  One  of  them  is  removing  the  general  position 
assumptions  and  defining  a  cell  of  a  set  of  sites  to  contain  all  points  of  P  to 
which  all  sites  of  the  set  are  simultaneously  closest  (see,  for  example,  [ES]). 
The  objective,  once  again,  would  be  to  estimate  the  complexity  of  this 
diagram  and  devise  an  efficient  algorithm  for  constructing  it.  Another  possi- 
ble avenue  for  extending  our  results  is  removing  the  restriction  that  sites  be 
points.  Can  our  approach  be  extended  to  yield  an  efficient  algorithm  for  com- 
puting the  diagram  of  a  set  of  straight-line  segments  in  a  polygon  under  a 
suitable  metric,  e.g.,  using  an  appropriate  modification  of  Yap's  technique 
[Ya]?  It  might  also  be  interesting  to  consider  which  properties  of  the  geodesic 
distance  (such  as  those  discussed  in  Section  3)  survive  if  we  allow  the  sides  of 
P  to  be  curvilinear. 

A  somewhat  unrelated  question  that  can  be  raised  in  conjunction  with  a 
geodesic  Voronoi  diagram  is  whether  the  sweeping  approach  of  Fortune  [F] 
can  be  generalized  for  constructing  such  diagrams.  A  significant  obstacle  to 
such  a  generalization  is  the  absence  of  an  obvious  natural  equivalent  of  a 
straight  sweep-line. 

The  natural  problem  of  analyzing  the  furthest-point  geodesic  Voronoi 
diagram  is  considered  in  [AFW].  They  show  a  &(n+k)  bound  on  the  com- 
plexity of  the  diagram  and  produce  an  0{(n  -l-^)log(n  -1-^))  algorithm  for  con- 
structing it. 
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